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The  degradation  and  ultimately  lifetime  of  solid  oxide  fuel  cells  (SOFCs)  is  determined  in  part  by  the 
stresses  generated  within  the  different  layers  of  the  device.  For  fully  dense  materials  such  as  the  elec¬ 
trolyte,  when  modelling  these  stresses  on  a  macro-scale  the  material  properties  can  be  considered  to  be 
homogeneous  (evenly  distributed)  allowing  the  prediction  of  volume  average  stresses  due  to  differential 
thermal  expansion  in  the  layer.  However,  detailed  stress  analysis  of  real,  multiphase  porous  layers  such 
as  those  found  in  SOFC  electrodes,  on  the  micron  and  sub-micron  scale  has  not  been  possible  to  date  as 
detailed  geometry  and  convenient  methods  to  generate  a  finite  element  model  have  not  been  available. 

In  this  paper  we  present  work  that  combines  microstructural  characterisation  of  a  porous  solid  oxide 
fuel  cell  anode  with  three  dimensional  stress  analysis  to  inspect  the  stresses  within  the  individual  phases 
of  the  anode,  and  at  phase  boundaries.  The  electrode  microstructure  has  been  characterised  using  focused 
ion  beam  (FIB)  tomography  and  the  resulting  microstructure  used  to  generate  a  solid  mesh  of  three 
dimensional  tetrahedral  elements.  A  temperature  field  was  applied  to  simulate  the  heating  of  the  sample 
from  room  temperature  (298  K)  to  operating  temperature  (1073  K).  The  maximum  principal  stress  in  the 
nickel  phase  was  found  to  exceed  the  yield  strength,  while  the  minimum  principal  stress  in  the  yttria- 
stabilized  zirconia  (YSZ)  phase  was  found  to  exceed  the  characteristic  strength  of  that  volume  of  YSZ, 
indicating  that  the  probability  of  failure  of  the  YSZ  matrix  is  significant. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  a  high  temperature  fuel  cell 
variant  that  employ  oxide  ion  conducting  ceramic  electrolytes 
operating  in  the  temperature  range  of  500-1 000  °C.  These  devices 
show  excellent  potential  for  a  broad  range  of  power  generation 
applications  and  particularly  for  large  scale  and  residential  com¬ 
bined  heat  and  power  (CHP)  [1].  However,  elevated  temperature 
operation  imparts  certain  challenges,  not  least  the  management  of 
mechanical  stresses  due  to  differential  thermal  expansion  between 
layers,  and  between  different  phases  within  the  electrode  layers. 

The  Ni-YSZ  (yttria-stabilized  zirconia)  cermet  is  a  common 
choice  for  the  anode  material  of  SOFCs  in  both  electrolyte  sup¬ 
ported  and  anode  supported  geometries  owing  to  its  low-cost, 
thermal  and  chemical  stability,  good  electronic  conductivity,  and 
high  catalytic  activity  for  hydrogen  oxidation  [2].  Providing  inti¬ 
mate  contact  of  ionic,  electronic  and  pore  phases  at  triple  phase 


*  Corresponding  author.  Tel.:  +44  0  20  7594  1539;  fax:  +44  0  20  7594  6566. 
E-mail  address:  r.clague@imperial.ac.uk  (R.  Clague). 

0378-7753 /$  -  see  front  matter  ©  201 1  Elsevier  B.V.  All  rights  reserved, 
doi:  1 0.1 01 6/j.jpowsour.201 0.1 2.097 


boundaries  (TPBs)  is  critical  to  optimising  electrochemical  perfor¬ 
mance.  In  these  composite  porous  anodes,  TPBs  are  distributed 
throughout  the  electrode  bulk. 

Recently,  the  use  of  advanced  focused  ion  beam  (FIB)  and  X- 
ray  nanofocus  computed  tomography  (nano-CT)  techniques  has 
enabled  direct  access  to  information  on  the  spatial  distribution  of 
TPBs  in  an  SOFC  electrode  structure  (e.g.  [3-10]).  These  tomographs 
provide  unprecedented  detail  regarding  the  complex  microstruc¬ 
ture  of  SOFC  electrodes  and  can  be  used  to  understand  evolution 
and  degradation  mechanisms  during  processing  and  operation. 
However,  these  tomographic  techniques  are  perhaps  most  pow¬ 
erful  when  utilised  in  tandem  with  numerical  simulation  to 
explicitly  evaluate  the  influence  of  microstructural  heterogeneities, 
as  demonstrated  previously  in  other  materials  ranging  from  metal 
matrix  composites  [11]  to  glass  foams  [12]. 

The  literature  provides  a  wealth  of  electrode  modelling  studies 
exploring  processing,  electrochemical  operation  and  mechanical 
performance.  These  models  typically  adopt  some  microstructural 
framework  of  varying  complexity  in  two  or  three  dimensions; 
indeed  there  has  been  significant  effort  in  developing  tools  for  gen¬ 
eration  of  microstructural  analogues  [13-15].  While  these  tools 
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are  increasingly  capable  of  generating  realistic  microstructural 
analogues  they  remain  idealised  representations.  The  ability  to 
interface  tomographic  techniques  with  relevant  simulations  is 
therefore  highly  desirable. 

A  number  of  groups  have  recently  demonstrated  this  capability: 
Shearing  et  al.  utilised  a  volume  of  fluid  (VOF)  model  to  predict  elec¬ 
trochemical  performance  in  a  Ni-YSZ  anode  [16].  Using  voxelised 
FIB  data,  the  microstructure  was  discretised  into  ca.  8000  elements, 
transport  phenomena  and  electrochemical  behaviour  was  mod¬ 
elled  in  each  of  the  VOF  elements  providing  a  prediction  of  the 
spatial  distribution  of  current  generation. 

Lattice  Boltzmann  methods  have  also  been  successfully  applied 
to  SOFC  reconstructions  from  FIB  tomography  [4]  and  X-ray  nano- 
CT  [7,9].  These  methods  enable  the  direct  use  of  voxelised  data 
without  the  requirement  for  “coarse-meshing”.  Flowever,  the  tech¬ 
niques  are  computationally  expensive,  often  requiring  access  to 
high  performance  super  computers  to  model  representative  vol¬ 
ume  elements.  Recently  Chiu  et  al.  estimated  that  using  a  linearly 
scaled  Lattice  Boltzmann  model  for  methane  reforming  in  a  real- 
life  microstructure  obtained  from  tomography  experiments  would 
take  up  to  8  weeks  on  a  2000  processor  facility  [17]. 

Wilson  et  al.  demonstrated  the  potential  for  combining  FE  mod¬ 
els  with  tomography  data  for  prediction  of  gas  phase  tortuosity 
[3].  Flowever,  the  community  has  typically  been  constrained  by  the 
ability  to  generate  high  quality  meshes  for  conventional  FE  analysis 
-  in  the  authors’  experience  meshing  tools  can  be  liable  to  artefacts 
(particularly  at  interfaces)  which  require  laborious  manual  repair. 
The  ability  to  accurately  mesh  microstructural  tomography  data  for 
combination  with  FE  simulations  will  provide  a  powerful  tool. 

Flere  we  present  the  combination  of  FIB  tomography  data  with 
mechanical  stress  simulation.  A  tetrahedral  finite  element  mesh 
was  successfully  created  from  the  reconstructed  tomography  sur¬ 
faces  using  commercially  available  software  (Scan  FE,  Simpleware 
Ltd.,  Exeter). 

2.  Microstructure  reconstruction 

2.1.  FIB  technique 

Electrolyte  supported  Ni-YSZ  symmetrical  cells  were  manu¬ 
factured  by  screen  printing  commercially  available  NiO-YSZ  ink 
(Ceramtec,  Germany)  onto  commercially  available  tape-cast  elec¬ 
trolyte  supports  (Nextec  Inc.  USA).  The  half  cells  were  subsequently 
chemically  reduced  and  characterised  using  electrochemical 
impedance  spectroscopy  ( see  [16])  Following  electrochemical  test¬ 
ing,  FIB  tomography  was  used  to  characterise  the  Ni-YSZ  electrode 
microstructure  using  a  Zeiss  xB-1540  FIB-FEG  SEM  at  the  London 
Centre  for  Nanotechnology. 

Using  the  dual  beam  FIB  to  sequentially  mill  and  image  the  Ni- 
YSZ  electrode  structure,  a  sequence  of  170  images  was  collected; 
pixel  dimensions  of  20  nm  were  recorded  with  an  inter-slice  thick¬ 
ness  of  15  nm.  From  this  sequence,  100  images  were  selected  and 
segmented  for  further  analysis.  Further  details  of  the  tomography 
data  collection  have  been  published  elsewhere  [16]. 

The  total  reconstructed  volume  is  6.68  x  5.04  x  1.50  p,m.  While 
larger  reconstructions  have  been  conducted  (e.g.  [5,8])  this  vol¬ 
ume  is  representative  of  the  screen-printed  electrode,  which  has  a 
thickness  of  5  p,m.  The  total  TPB  density  measured  following  a  tech¬ 
nique  described  by  Shearing  et  al.  [16]  is  12.99  |jnrr2  with  volume 
fractions  Ni  =  35.5%,  YSZ  =  46.7%  and  pore  =  1 7.8%. 

3.  Finite  element  model 

The  finite  element  (FE)  model  was  generated  from  the  seg¬ 
mented  tomography  data  using  ScanFE  software  (ScanFE  V3.1, 


YSZ  Phase  mesh 


Ni  Phase  mesh 


Fig.  1.  (a)  YSZ  phase  mesh,  (b)  Ni  phase  mesh 


Simpleware  Ltd,  UK).  The  FE  model  consisted  of  285,000  four-node 
tetrahedral  elements  (Abaqus  type  C3D4).  The  appropriate  num¬ 
ber  of  elements  was  determined  from  a  mesh  independence  test, 
which  showed  that  mesh  independent  stress  results  based  on  aver¬ 
age  stress  were  obtained  at  this  element  density.  During  meshing 
the  Ni  and  YSZ  phases  were  assigned  to  one  of  two  material  sections 
to  allow  allocation  of  appropriate  material  properties,  as  shown  in 
Fig.  1(a)  and  (b).  Node  and  element  sets  were  automatically  gen¬ 
erated  on  the  model  cut  planes  and  at  the  interface  of  the  two 
phases. 

Boundary  conditions  and  load-cases  in  the  form  of  a  global  tem¬ 
perature  field  were  applied  (as  detailed  below)  using  Abaqus  CAE 
(v6.10,  Simulia,  UK)  and  finally  a  stress  analysis  was  run  using  the 
Abaqus  Implicit  finite  element  solver.  The  analysis  assumed  that 
there  was  no  separation  of  the  Ni  and  YSZ  phases  and  that  the 
material  behaviour  was  linear  and  elastic. 

The  simulation  took  89  s  to  complete  on  a  computer  with  a  Linux 
operating  system  using  four  AMD  Opteron  dual  core  64  bit  proces¬ 
sors  and  a  2  GHz  clock  speed. 

3.1.  Boundary  conditions 

In  order  to  represent  the  reconstructed  portion  of  the  electrode 
as  a  section  of  a  larger  anode,  the  boundary  conditions  that  were 
applied  to  the  model  constrained  the  cut  planes  to  remain  pla- 
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Table  1 

Material  properties. 


Layer 

Material 

Temperature  (K) 

E 

V 

a(x  10-6) 

Yield  (MPa) 

UTS  (MPa) 

Anode 

Ni  [18] 

1073 

207 

0.31 

13.5 

59 

317 

YSZ  [19] 

1073 

157 

0.313 

10.5 

n/a 

n/a 

nar,  while  allowing  bulk  thermal  expansion  between  opposite  cut 
planes.  This  was  achieved  by  restraining  the  degree  of  freedom  of 
the  nodes  on  one  cut  plane  in  the  direction  normal  to  the  cut  plane, 
then  on  the  opposite  cut  plane  the  normal  degree  of  freedom  of 
the  nodes  on  that  plane  were  tied  together.  This  has  the  effect 
of  allowing  bulk  thermal  expansion  between  opposing  cut  planes, 
while  still  imposing  the  condition  that  the  cut  planes  remain  pla¬ 
nar.  These  mechanical  boundary  conditions  are  not  considered  to 
fully  represent  the  behaviour  of  the  cut  planes  of  the  sample  were 
it  in  a  fuel  cell,  but  this  will  be  investigated  in  due  course  through 
sub-modelling. 

As  the  YSZ  phase  is  wetted  by  the  Ni  phase,  a  strain  continu¬ 
ity  condition  was  imposed  at  the  interface  between  the  Ni  and 
YSZ  phases.  The  solid  phases  were  free  to  expand  into  the  pore 
phase.  A  global  temperature  field  T0  of  298  K  was  applied  as  an  ini¬ 
tial,  zero  stress  condition,  which  was  subsequently  increased  to  a 
temperature  T\  of  1073  K  in  a  single  analysis  step. 

3.2.  Material  properties 


a 


The  mechanical  properties  of  the  materials  required  for  the 
thermo-mechanical  stress  analysis  are  summarised  in  Table  1.  E 
is  the  elastic  modulus  (GPa)  v  is  the  dimensionless  Poisson’s  ratio 
and  a  is  the  coefficient  of  thermal  expansion  (1  K-1 ). 


YSZ  Phase 
Ni  Phase 


Distribution  of  phases  in  sample 


4.  Results 

4.1.  Stress  analysis 

The  differential  thermal  expansion  of  the  Ni  and  YSZ  phases 
generates  stresses  within  the  composite  structure  that  peak  at  the 
interfaces  of  the  solid  phases.  Figs.  2(a),  (b)  and  3  shows  the  maxi¬ 
mum  principal  stress  field  generated  when  the  global  temperature 
of  the  composite  structure  changes  from  T0  to  T\.  In  the  Ni  phase 
the  peak  maximum  principal  stress  of  over  200  MPa  is  a  significant 
proportion  of  the  failure  strength  of  31 7  MPa,  and  well  above  the 
yield  strength  of  57  MPa  at  1 073  K,  illustrating  that  the  elastic  anal¬ 
ysis  used  is  a  simplifying  approximation  that  will  over-estimate 
peak  stresses. 

5.  Discussion 

The  stresses  predicted  by  this  analysis  exceed  the  yield  strength 
of  the  Ni  phase  by  some  margin,  implying  that  stress  relief  through 
plastic  deformation  is  likely  to  occur  at  the  phase  boundaries  as  the 
sample  is  heated  from  room  to  operating  temperature.  By  assuming 
perfectly  elastic  material  behaviour  and  a  zero  stress  state  at  room 
temperature,  the  analysis  predicts  that  plastic  deformation  would 
occur  within  the  Ni  phase  near  the  phase  boundary.  Unfortunately 
the  strength  of  the  Ni-YSZ  interface  is  not  known  so  no  comment 
can  be  made  on  whether  delamination  is  expected. 

In  the  YSZ  phase  the  minimum  principal  stress  has  a  magni¬ 
tude  of  over  200  MPa.  The  failure  strength  of  YSZ  is  not  a  simple 
matter  to  determine  as  it  is  a  brittle  material  and  therefore  the  fail¬ 
ure  strength  is  volume  dependent  [20].  It  is  more  usual  to  consider 
probability  of  failure  when  talking  about  brittle  materials  but  in 
this  analysis,  in  order  to  qualitatively  examine  whether  the  stresses 
are  more  likely  than  not  to  cause  failure  the  stresses  are  simplisti- 
cally  compared  to  the  characteristic  strength,  ae,  of  the  YSZ  phase, 


b 


Section  A-A 


Maximum  Principal  Stress 

Fig.  2.  (a)  Distribution  of  phases,  (b)  maximum  principal  stress 
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3D  Max  Principal  Stress 


Fig.  3.  3D  max  principal  stress 

defined  as  the  stress  at  which  50%  of  samples  of  this  structure  would 
fail.  For  the  purposes  of  this  analysis  the  volume  of  material,  Ve  was 
taken  to  be  equal  to  the  volume  of  the  reconstructed  sample  and 
ae  was  calculated  from  Equation  (1).  The  material  property  a0,  is 
defined  as  the  stress  at  which  samples  of  the  material  of  unit  volume 
would  fail.  The  values  used  here  are  taken  as  g0  =  23.1  MPa  [21  ],  the 
Weibull  modulus  m  =  8.6  [19],  and  Ve  =  50  x  10-9  m3,  giving  a  value 
of  Gq  =  163  MPa. 

G0  =  GO{Ve)~^  (1) 

Assuming  the  failure  characteristics  of  YSZ  are  the  same  in 
compression  and  tension,  the  characteristic  strength  of  YSZ  at 
these  volumes  is  therefore  exceeded,  so  from  the  results  of 
this  analysis  some  fracturing  of  the  ceramic  matrix  would  be 
expected.  In  this  analysis  linear  elastic  material  properties  have 
been  assumed,  which  will  tend  to  over  predict  stress  levels, 
as  plastic  deformation  in  the  nickel  phase  will  tend  to  relieve 
stress.  A  time  dependent  elastic-plastic  analysis  of  the  structure 
is  expected  to  predict  lower  final  stress  levels  and  therefore  a 
lower  probability  of  failure.  While  a  calculation  of  the  probabil¬ 
ity  of  failure  of  the  structure  is  beyond  the  scope  of  this  paper, 
using  a  method  demonstrated  by  Clague  [21]  the  probability  of 
failure  of  each  finite  element  of  the  ceramic  phase  can  be  calcu¬ 
lated,  as  the  mean  stress  and  volume  of  each  finite  element  are 
known.  The  probability  of  failure  for  the  entire  phase  can  then 
be  calculated  by  summing  the  individual  failure  probabilities  of 
the  finite  elements  following  the  Weibull  weakest  link  theory 
[20]. 

6.  Conclusion 

This  work  has  demonstrated  the  viability  of  a  technique  to  cap¬ 
ture  real,  highly  complex,  two  phase  microstructural  geometry 
such  as  that  found  in  SOFC  electrodes  using  FIB  tomography  and 
generate  a  finite  element  model  of  suitable  quality  for  thermo¬ 
mechanical  stress  analysis  from  it.  Using  simple  elastic  analysis, 


the  method  was  shown  to  be  a  potentially  powerful  tool  for  study¬ 
ing  the  response  of  complex  multiphase  structures  to  thermally, 
mechanically  or  chemically  induced  stress.  The  key  result  in  this 
work  is  that  the  yield  strength  of  nickel  is  exceeded  at  the  inter¬ 
face  of  the  Ni  and  YSZ  phases,  indicating  that  this  is  potentially  an 
important  mechanism  for  stress  relief  in  SOFC  electrodes  as  they 
are  heated  or  cooled. 

The  authors  are  presently  developing  a  technique  to  examine 
the  normal  stresses  at  phase  boundaries,  allowing  a  prediction  of 
TPB  movement  or  loss  when  compared  to  the  work  of  adhesion  of  Ni 
on  YSZ.  This  technique  offers  the  possibility  of  studying  the  stresses 
within,  and  failure  mechanisms  of,  the  nickel  and  YSZ  phases.  This 
method  also  has  application  as  a  submodelling  technique,  with  dis¬ 
placement  boundary  conditions  applied  at  the  model  cut  planes 
from  a  macro-scale  stress  analysis. 
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